**************************************************************
* Monte-Carlo: LP-DiD vs XTEVENT  (event time −3 … +3)
**************************************************************
clear all
* RUN 00_path_master.do FIRST TO GET FILEPATHS
set seed 60112

*--------------------------- 1. globals ----------------------------
local simulations   = 500
local leads         = 3
local lags          = 3
local total_coefs   = `leads' + `lags' + 1      // 7 columns: m3…p3

local states        = 50
local time_periods  = 15
local never_treated = 20
local early_treated = 12
local late_treated  = 18

cap frame drop sim*

*-------------------- 2. Monte-Carlo (simtype 1–4) -----------------
forval simtype = 4/4 {

    *---- storage matrices (2 estimators × 7 coefs) ----------------
    matrix results_lpdid   = J(`simulations',`total_coefs',.)
    matrix results_xtevent = J(`simulations',`total_coefs',.)
    matrix t_lpdid         = J(`simulations',`total_coefs',.)
    matrix t_xtevent       = J(`simulations',`total_coefs',.)

    *------------------ 2A. simulation loop ------------------------
    forval sim = 1/`simulations' {
        di "Sim `sim'"
        /* ----------  data-generating process -------------------- */
        clear
        set obs `states'
        gen state_id = _n
        gen group    = cond(state_id<=`never_treated',0, ///
                       cond(state_id<=`never_treated'+`early_treated',1,2))
        expand `time_periods'
        bys state_id: gen time = _n-1

		* Draw one state FE per state, then carry it down each state's panel
		bysort state_id (time): gen stateFE = cond(_n==1, rnormal(0,.15), .)
		bysort state_id (time): replace stateFE = stateFE[1]
		
		* Draw one time FE per calendar year, then carry it across states
		bysort time (state_id): gen timeFE = cond(_n==1, rnormal(0,.15), .)
		bysort time (state_id): replace timeFE = timeFE[1]
        gen trend = 0
        if `simtype'==1 replace trend = -.25*time if inlist(group,1,2)
        if `simtype'==1 replace trend = -.25*5    if inlist(group,1,2) & time>=5
        if `simtype'==2 replace trend = -.25*time if inlist(group,1,2)
        if `simtype'==2 replace trend = -.25*5    if inlist(group,1,2) & time>=5
        if `simtype'==3 replace trend = -.25*time if group==1
        if `simtype'==3 replace trend = -.25*5    if group==1 & time>=5
        if `simtype'==4 replace trend = -.2*time  if group==1
        if `simtype'==4 replace trend =  .9*time  if group==2
        if `simtype'==4 replace trend = -.2*5     if group==1 & time>=5
        if `simtype'==4 replace trend =  .9*5     if group==2 & time>=5
        replace time = time + 1

        gen treatment = 0
        replace treatment = 1 if (group==1 & time>=9)
        replace treatment = 1 if (group==2 & time>=12)
        if `simtype'==1 {
            replace treatment = 0 if group==2
            replace treatment = 1 if (group==2 & time>=9)
        }

        gen cohort  = cond(group==1,9, cond(group==2,12,0))
        if `simtype'==1 replace cohort = 9 if group==2
        xtset state_id time

        *-------- lead/lag dummies (≤−3 binned) --------------------
        forval j = 1/2 {
            gen lead`j' = D.F`j'.treatment   // F1, F2
        }
        gen lead3plus = F3.treatment
        forval j = 4/15 {
            cap replace lead3plus = 1 if F`j'.treatment==1
        }

        forval j = 1/3 {
            gen lag`j' = D.L`j'.treatment    // L1, L2, L3
        }
        gen dtreatment = D.treatment
        replace dtreatment = 0 if dtreatment==.

        gen error   = rnormal(0,.15)
        gen outcome = trend + error + stateFE + timeFE

        /******************** LP-DiD (–3+ bin) *********************/
        gen y_m1 = L1.outcome
        gen y_m2 = L2.outcome
        gen y_p0 = outcome
        gen y_p1 = F1.outcome
        gen y_p2 = F2.outcome
        gen y_p3 = F3.outcome

        bys state_id (time): gen csum_y = sum(outcome)
        bys state_id (time): gen csum_n = _n
        gen y_m3plus = L3.csum_y / L3.csum_n      // mean k≤−3
		
		
		bys state_id       : gen  total_y = csum_y[_N]       // state-total
		bys state_id (time): gen  rev_csum_y = total_y - (csum_y - outcome)  // t→end
		bys state_id (time): gen  rev_csum_n = _N - _n + 1                 // obs count
		gen y_p3plus = F3.rev_csum_y / F3.rev_csum_n    // k ≥ +3   (fixed)
		gen y_p2plus = F2.rev_csum_y / F2.rev_csum_n    // k ≥ +3   (fixed)


        gen dif_m3 = y_m3plus - y_m1
        gen dif_m2 = y_m2     - y_m1
        gen dif_m1 = 0
        gen dif_p0 = y_p0     - y_m1
        gen dif_p1 = y_p1     - y_m1
        gen dif_p2 = y_p2     - y_m1
        gen dif_p3 = y_p3plus - y_m1

        gen sample_lp = dtreatment==1 | group==0   // newly treated + never

        local diflist "dif_m3 dif_m2 dif_p0 dif_p1 dif_p2 dif_p3"
        local collist  1       2       4       5       6       7
        local i = 1
        foreach d of local diflist {
            local c = word("`collist'",`i')
            qui reghdfe `d' dtreatment if sample_lp, absorb(time) vce(cluster state_id)
            matrix results_lpdid[`sim',`c'] = _b[dtreatment]
            matrix t_lpdid[`sim',`c']       = abs(_b[dtreatment]/_se[dtreatment])
            local ++i
        }
        matrix results_lpdid[`sim',3] = 0       // k = –1
        matrix t_lpdid[`sim',3]       = 0

        /******************** XTEVENT (w = 2) **********************/
        xtevent outcome, w(2) pol(treatment) nofe ///
                 cluster(state_id) reghdfe addabsorb(state_id time)

        matrix results_xtevent[`sim',1] = _b[_k_eq_m3]
        matrix results_xtevent[`sim',2] = _b[_k_eq_m2]
        matrix results_xtevent[`sim',3] = 0
        matrix results_xtevent[`sim',4] = _b[_k_eq_p0]
        matrix results_xtevent[`sim',5] = _b[_k_eq_p1]
        matrix results_xtevent[`sim',6] = _b[_k_eq_p2]
        matrix results_xtevent[`sim',7] = _b[_k_eq_p3]

        matrix t_xtevent[`sim',1] = abs(_b[_k_eq_m3]/_se[_k_eq_m3])
        matrix t_xtevent[`sim',2] = abs(_b[_k_eq_m2]/_se[_k_eq_m2])
        matrix t_xtevent[`sim',3] = 0
        matrix t_xtevent[`sim',4] = abs(_b[_k_eq_p0]/_se[_k_eq_p0])
        matrix t_xtevent[`sim',5] = abs(_b[_k_eq_p1]/_se[_k_eq_p1])
        matrix t_xtevent[`sim',6] = abs(_b[_k_eq_p2]/_se[_k_eq_p2])
        matrix t_xtevent[`sim',7] = abs(_b[_k_eq_p3]/_se[_k_eq_p3])
    }  /* ---- end simulation loop ------------------------------ */

    *------------------ 3. matrix → variables ---------------------
    svmat results_lpdid,   names(lp)
    svmat results_xtevent, names(xtevent)
    svmat t_lpdid,   names(t_lp)
    svmat t_xtevent, names(t_xtevent)

    gen lag = _n-4     // -3 -2 -1 0 1 2 3

    /* -------- 4. share |t| ≥ 1.96 ----------------------------- */
    local estlist "t_lp t_xtevent"
    matrix share = J(2,`total_coefs',.)
    local r = 1
    foreach est of local estlist {
        forval j = 1/`total_coefs' {
            gen _sig = (abs(`est'`j') >= 1.96)
            summarize _sig, meanonly
            matrix share[`r',`j'] = r(mean)
            drop _sig
        }
        local ++r
    }
    matrix rownames share = LPDID XTEVENT
    matrix colnames share = m3 m2 m1 p0 p1 p2 p3
    mat list share, format(%5.3f)

    /* -------- 5. means & 95% CI (for plot) -------------------- */
    foreach pref in lp xtevent {
        gen coef_mean_`pref' = .
        gen coef_sd_`pref'   = .
        forval i = 1/7 {
            summarize `pref'`i', detail
            replace coef_mean_`pref' = r(mean) in `i'
            replace coef_sd_`pref'   = r(sd)   in `i'
        }
        gen `pref'_lower = coef_mean_`pref' - 1.96*coef_sd_`pref'
        gen `pref'_upper = coef_mean_`pref' + 1.96*coef_sd_`pref'
    }

}   /* ------------- end simtype loop ---------------------- */


**************************************************************
* Figures  (DGP & Simulation results)
**************************************************************
gen _outcome = outcome - error - stateFE - timeFE
replace _outcome = _outcome + 0.03 if state_id == 1
replace _outcome = _outcome - 0.03 if state_id == 41
replace _outcome = _outcome + 0.03 if state_id == 31

tw (line _outcome time if state_id==1,  lcol(blue%70)  lpat(dash)) ///
   (line _outcome time if state_id==31, lcol(red%70)   lpat(dash)) ///
   (line _outcome time if state_id==31 & treatment, lcol(red%40)   lpat(sol) lw(1.5)) ///
   (line _outcome time if state_id==41, lcol(green%70) lpat(dash)) ///
   (line _outcome time if state_id==41 & treatment, lcol(green%40) lpat(sol) lw(1.5)), ///
   ytitle("Mean Simulated Outcome") xtitle("Calendar time") ///
   legend(order(1 "Never-treated" 2 "Early adopters (no treat)" 3 "Early adopters (treated)" ///
                4 "Late adopters (no treat)" 5 "Late adopters (treated)") row(2) pos(6)) ///
   ysize(5) xsize(9)
graph export "$figures/Figure C1_A.pdf", replace

  
	   
twoway                                                          ///
    /* ─── LP-DiD ribbon & lines ───────────────────────────── */ ///
    (rarea  lp_lower  lp_upper  lag if inrange(lag,-3, 3),      ///
            fcolor(forest_green%20) lcolor(forest_green%0))      ///
    (line    coef_mean_lp lag if inrange(lag,-2, 3),            ///
            lcolor(forest_green%50) lwidth(med))                ///
    (line    coef_mean_lp lag if inrange(lag,-3,-2),            ///
            lcolor(forest_green%50) lwidth(med) lpattern(dash)) ///
    (scatter coef_mean_lp lag if lag==-3,                       ///
            msymbol(Th) mcolor(forest_green))                   ///
    (scatter coef_mean_lp lag if inrange(lag,-2, 3),            ///
            msymbol(triangle) mcolor(forest_green))             ///
    /* ─── XTEVENT ribbon & lines ──────────────────────────── */ ///
    (rarea  xtevent_lower xtevent_upper lag if inrange(lag,-3,3), ///
            fcolor(dkorange%20) lcolor(dkorange%0))             ///
    (line    coef_mean_xtevent lag if inrange(lag,-2, 3),       ///
            lcolor(dkorange%50) lwidth(med))                    ///
    (line    coef_mean_xtevent lag if inrange(lag,-3,-2),       ///
            lcolor(dkorange%50) lwidth(med) lpattern(dash))     ///
    (scatter coef_mean_xtevent lag if lag==-3,                  ///
            msymbol(Dh) mcolor(dkorange))                      ///
    (scatter coef_mean_xtevent lag if inrange(lag,-2, 3),       ///
            msymbol(D)  mcolor(dkorange))                      ///
    /* ───  brace (three short line segments)  ─────────────── */ ///
    (scatteri -0.65 -3.5   -0.65 -2.2,      recast(line) lwidth(vthin) lpat(solid) lcolor(forest_green)) ///
    (scatteri -0.65 -3.5   -0.6 -3.5,      recast(line) lwidth(vthin)  lpat(solid) lcolor(forest_green)) ///
    (scatteri -0.65 -2.2   -0.6 -2.2,      recast(line) lwidth(vthin) lpat(solid)  lcolor(forest_green)) ///
	/* ───  brace (three short line segments)  ─────────────── */ ///
    (scatteri -0.4 1.2   -0.4 3.2,      recast(line) lwidth(vthin) lpat(solid) lcolor(forest_green)) ///
    (scatteri -0.4 1.2   -0.35 1.2,      recast(line) lwidth(vthin)  lpat(solid) lcolor(forest_green)) ///
    (scatteri -0.4 3.2   -0.35 3.2,      recast(line) lwidth(vthin) lpat(solid)  lcolor(forest_green)), ///
    /* ───  label below brace  ─────────────────────────────── */ ///
    text(-0.65  -2.8  "LP-DiD correctly detects" "early period parallel-" "trends violations," ///
            , size(5pt) place(s) justification(center) color(forest_green))       ///
	text(-0.76  -2.7  " XTEVENT doesn't detect them" ///
            , size(5pt) place(s) justification(center) color(dkorange))  ///
    /* ───  label below brace  ─────────────────────────────── */ ///
    text(-0.42  2.2  "LP-DiD correctly detects zero effect," ///
            , size(5pt) place(s) justification(center) color(forest_green))       ///
	text(-0.47  2.2  "XTEVENT mistakenly detects" " negative effect" ///
            , size(5pt) place(s) justification(center) color(dkorange))  ///
    /* ───  axes, legend, sizing  ──────────────────────────── */ ///
      xtitle("Event date")                                     ///
       xlab(-3 "{&le}-3" -2 "-2" -1 "-1" 0 "0" 1 "1" 2 "2"      ///
            3 "{&ge}3")                                         ///
       ytitle("Cumulative response")                            ///
       legend(order(10 "XTEVENT (TWFE)" 5 "LP-DiD") row(1) pos(6))     ///
       ysize(3) xsize(4)

graph export "$figures/Figure C1_B.pdf", replace


******************************************************************
* 7.  Summary matrices & LaTeX table  (pretty version)
******************************************************************
local est_prefix "lp xtevent"            // 1 = LP-DiD, 2 = XTEVENT
local est_names  "LPDID XTEVENT"
local total_coefs = 7                    // m3 … p3

matrix pos    = J(2,`total_coefs',.)
matrix neg    = J(2,`total_coefs',.)

local R = 0
foreach e of local est_prefix {
    local ++R
    forvalues j = 1/`total_coefs' {
        count if !missing(`e'`j')
        local ntot = r(N)
        count if abs(t_`e'`j')>=1.96 & `e'`j' > 0 & `e'`j' <.
        matrix pos[`R',`j'] = r(N)/`ntot'
        count if abs(t_`e'`j')>=1.96 & `e'`j' < 0
        matrix neg[`R',`j'] = r(N)/`ntot'
    }
}
matrix rownames pos = `est_names'
matrix rownames neg = `est_names'

matrix colnames pos  = m3 m2 m1 p0 p1 p2 p3
matrix colnames neg  = m3 m2 m1 p0 p1 p2 p3

*------------------------------------------------------------------
* 8.  Write LaTeX table
*------------------------------------------------------------------
tempname fh
file open `fh' using "$tables/Table C1.tex", write replace text
file write `fh' "{\small" _n
file write `fh' "\begin{tabular}{l*{7}{c}}" _n
file write `fh' "\toprule" _n
file write `fh' " & \multicolumn{7}{c}{Coefficients for event times}\\\cmidrule(lr){2-8}" _n
file write `fh' " & \$\leq-3\$ & \$-2\$ & \$-1\$ & \$0\$ & \$1\$ & \$2\$ & \$\geq3\$\\\midrule" _n

* helper to print a row
cap program drop _rowout
program define _rowout
        args MAT ROW LABEL HANDLE
        local line "`LABEL'"
        forvalues j = 1/7 {
                local cell : display %6.3f `MAT'[`ROW',`j']
                local line "`line' & `cell'"
        }
        file write `HANDLE' "`line' \\" _n
end

* -------- Panel A : LP-DiD --------------------------------------
file write `fh' "\multicolumn{8}{l}{\textbf{Panel A: LP-DiD}}\\\\\addlinespace" _n
_rowout pos 1 "Positive and sig." `fh'
_rowout neg 1 "Negative and sig." `fh'
file write `fh' "\addlinespace" _n

* -------- Panel B : XTEVENT ------------------------------------
file write `fh' "\multicolumn{8}{l}{\textbf{Panel B: XTEVENT}}\\\\\addlinespace" _n
_rowout pos 2 "Positive and sig." `fh'
_rowout neg 2 "Negative and sig." `fh'

file write `fh' "\bottomrule" _n
file write `fh' "\end{tabular}}" _n
file close `fh'

di as text "New LaTeX table written to Table C1.tex"

